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Some of the most popular ways to treat quantum critical materials, that is, materials close to a 
magnetic instability, are based on the Landau functional. The central quantity of such approaches 
is the average magnitude of spin fluctuations, which is very difficult to measure experimentally or 
compute directly from the first principles. We calculate the parameters of the Landau functional for 
Pd and use these to connect the critical fluctuations beyond the local-density approximation and 
the band structure. 



The physics and materials science of weak itinerant fer- 
romagnetic metals and highly renormalized paramagnets 
near magnetic instabilities has attracted renewed theo- 
retical interest. This is a result of recent discoveries of 
materials with highly non-conventional metallic proper- 
ties, especially, non-Fermi liquid scalings, metamagnetic 
behavior, and unconventional superconductivity, in sev- 
eral cases co-existing with ferromagnetism. Discoveries 
in the last three years alone include the co-existing fer- 
romagnetism and superconductivity of ZrZn2 [1], UGe2 
[2],URhGe 2 [3], high pressure e-Fe [5], and the metamag- 
netic quantum critical point in Si'3Ru207 [4]. 

Unfortunately, although model theories have been put 
forth, there is still not an established material specific 
(first principles) theoretical understanding of these phe- 
nomena. One difficulty is the usual starting point for first 
principles theories, density functional theory (DFT) as 
implemented in the local density approximation (LDA). 
This already includes most spin degrees of freedom, in- 
cluding dynamical fluctuations, as evidenced by its for- 
mally exact description of the uniform electron gas as well 
as its well documented success in accurately describing 
a wide variety of itinerant magnetic materials. However, 
the electron gas, upon which most density functionals 
are built, is not near any critical point for densities rele- 
vant to the solid state, and furthermore the proximity to 
itinerant magnetism of a metal is an extremely non-local 
quantity, in particular depending on the electronic den- 
sity of states at the Fermi level N(Ep). Therefore, the 
exact DFT, which by definition includes all fluctuations 
and describes the ground state magnetization exactly, is 
likely to be extremely nonlocal and probably nonanalyt- 
ical for the materilas near a quantum critical point. 

On the other hand, the LDA, while providing a good 
description of most itinerant ferromagnets that are not 
near critical points, fails to include the soft critical fluctu- 
ations in the materials of interest here. Since fluctuations 
are generically antagonistic to ordering, the result is that 
magnetic moments and magnetic energies of weak itiner- 
ant ferromagnets near critical points are overestimated in 
the LDA, as opposed to LDA's failure to describe Mott- 
Hubbard insulators where the LDA underestimates the 
tendency to magnetism. Recent examples include Scsln 



[6], Ni 3 Al [13], NaCo 2 4 [8], and ZrZn 2 [9]. Similarly, 
susceptibilities of paramagnets near critical points are 
underestimated. Furthermore, there is an overlap re- 
gion where the LDA predicts ferromagnetism for para- 
magnetic materials. This interesting class includes FeAl 
[11,12], Ni 3 Ga [13], and Sr 3 Ru 2 7 [7] (as mentioned, this 
latter material shows a metamagnetic quantum critical 
point). The basic theoretical difficulty in correcting the 
LDA for these materials is that there is some unknown 
and possibly strongly material dependent cross-over in 
energy (and possibly non-trivially in momentum) sepa- 
rating quantum critical fluctuations, not included in the 
LDA, from the dynamical fluctuations that are included 
in the LDA. Qualitatively, this may be understood from 
the fact that the LDA is based on the properties of the 
uniform electron gas, which is far from any magnetic crit- 
ical point at densities relavent for solids, the consequence 
being a mean-field-like description of magnetism near 
critical points. Thus the underlying reason for the fail- 
ure of the LDA to describe these systems is very different 
from the failures in the well-known class of Coulomb cor- 
related materials, such as the Mott-Hubbard insulators. 
There, the basic problem is the neglect of some electron- 
electron interactions, and can often be largely corrected 
at the static level, e.g. via approaches like LDA+U. It 
is worth noting that these dynamical fluctuations are re- 
sponsible not only for the suppression of the magnetic 
ordering, but also for unusual transport properties of 
quantum critical materials, deviating from the conven- 
tional Fermi liquid behavior, for mass renormalization, 
and even for superconductivity in some systems. Many 
of these issues have been addressed recently in theoret- 
ical papers, utilizing idealized models of various kinds. 
However, a quantitative link between such models and 
actual material characteristics is still missing. 

We attempt to build a bridge between such theo- 
ries and the LDA. We concentrate on the question of 
what kind of material-specific understanding, relevant 
for quantum criticality, can be extracted from the LDA 
calculations. Primarily, we focus here on Pd. This is 
perhaps the best studied high susceptibility paramagnet 
[14-17], and in fact a number of theories related to spin 
fluctuations have been elucidated using this material. 
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Furthermore, itinerant ferromagnetism appears in Pd at 
2.5% Ni doping [18]. We present highly accurate calcu- 
lations of the static magnetic susceptibility for Pd and 
find that, indeed, the LDA overestimates the tendency 
to magnetism. We also estimate the r.m.s. magnitude of 
spin fluctuations (paramagnons) in Pd, needed to reduce 
the calculated susceptibility to reproduce experiment, 
and show that it is compatible with that which might be 
estimated from LDA susceptibility via the fluctuation- 
dissipation theorem with a reasonable ansatz for the cut- 
off momentum. 

We have performed electronic structure calculations 
using the self consistent full potential linearized aug- 
mented plane wave (FLAPW) [19] method within the 
density functional theory (DFT) [20]. The local density 
approximation (LDA) of Pcrdew and Wang [21] and the 
Generalized Gradient Approximation (GGA) of Perdew, 
Burke, and Ernzerhof [22] were used for the correlation 
and exchange potentials. Calculations were performed 
using the WIEN2k package [23] . Local orbital extensions 
[24] were included in order to accurately treat the upper 
core states and to relax any residual linearization errors. 
A well converged basis consisting of LAPW basis func- 
tions with wave vectors up to if max set as RK ma ^ = 9, 
with the Pd sphere radii R =2.59 bohr. All total energy 
calculations used at least 1470 and up to 2844 k-points 
in the irreducible part of the Brillouin zone as needed. 
Spin-orbit (SO) interactions were incorporated using a 
second variational procedure [25] , where all states below 
the cutoff energy 1.5 Ry were included, with the so-called 
Pi/2 extension [19], which accounts for the finite charac- 
ter of the wave function at the nucleus for the P1/2 state. 

All calculations were performed in an external mag- 
netic field, interacting with both spin, s, and orbital, 1, 
[23] momenta: 

V Hext = M B H ext (l + 2s). 

The input values of H were chosen from to 10000 T in 
irregular increments to map out the change in energy and 
magnetic moment as a function of applied field. While 
use of the LDA [21] resulted in zero magnetic moment 
in a zero magnetic field, consistent with the experiment 
[18], use of GGA [22] resulted in a persistent magnetic 
moment of 0.2/is, with an extremely small magnetic en- 
ergy of less than 1 meV. 

In order to understand the change in the total energy 
and magnetic moments as a function of the applied ex- 
ternal field, special care was taken to ensure that these 
quantities were well converged with respect to the k- 
mesh. Given that in the low fields we are interested in, 
energy changes need to be converged of the order of 0.1 
meV/atom. The total energy, E, with respect to that 
at M = hb as a function of the magnetization, M, is 
shown in Figure 1. Figure 2 shows the applied magnetic 
field, H, as a function of M (with the magnetization di- 
rection 100). Note that the latter dependence follows 



from the former one, as H = j^. One can see though 
that of the two quantities H shows less computational 
noise, so this was the dependency we used in the analysis 
described below. 

As can be seen in both plots (more so in Figure 2), 
there exist two regimes in terms of the magnetic moment, 
M. For values of M < 0.5 \ib (corresponding to H <~ 
1200 T), the external field and energy increase slowly, but 
for M > 0.5^b, both H and E increase rapidly, suggesting 
that the long wave spin fluctuations at any temperature 
should be smaller that ~ 0.5 /is in amplitude. 

The linear magnetic susceptibility is defined as = 
§^j\m=o — Figure 3 shows, however, that even for 

M < 0.5 hb the susceptibility is highly nonlinear. In 
fact, §^ starts near 11.6xl0~ 4 emu/mol and decreases 
rapidly with the field. In order to compute accurately the 
relevant derivatives, we have fitted the calculated H(M) 
for M < 0.5 hb with a polynomial (Figure 3). Thus com- 
puted susceptibility as the function of the applied field is 
shown in Fig. 4. We see that the zero field susceptibility 
is nearly twice larger than the experimental value of 6.8 
xl0~ 4 emu/mol corresponding to 21 st/eV-cell [30,31]. 
Only in a field of 550 T does the susceptibilty eventually 
become close to the experimental number. 

One may understand the origin of this overestimation 
of magnetic susceptibility in the following way. Not only 
is the calculated susceptibility very large, but also as 
mentioned the dependence of the induced magnetic mo- 
ment on the applied field is highly nonlinear in such a 
manner that the total energy as a function of the con- 
strained magnetic moment is very flat up to M « 0.5 [Ib- 
This implies that zero temperature quantum fluctuations 
beyond the LDA may have a substantial magnitude. One 
of the ways to take into account these fluctuations is via 
the Ginzburg-Landau theory, which, in connection with 
the spin fluctuations in nearly-magnetic metals has been 
used by several authors during the 1970's. This method 
starts with an expression for the total energy without 
such fluctuations as a function of the induced magnetic 
moment M 

E s taUc(M) = a + J2 ^ a2 « M2 "' W 

n>l 

Hstatic(M) — ^a 2n M 2n_1 (2) 

n>l 

(obviously, a2 gives the inverse spin susceptibility with- 
out fluctuations), and then assume Gaussian zero-point 
fluctuations of an r.m.s. magnitude £ for each of the d 
components of the magnetic moment (for a 3D isotropic 
material like Pd, d = 3). After averaging over the spin 
fluctuations, one obtains a fluctuation-corrected func- 
tional. The general expression of Ref. [32] can be written 
in the following compact form: 
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H(M) = ]T 5 2n M 2 "- 1 



2fc s 



a 2 „ = ECti^+ofitr 1 !! + T )- (3) 



i>0 



For instance, 



35 



35 



a 2 =a 2 + -a 4 ^ + y a e£ + y a s£ ••• 
a 4 = a 4 + y a 6^ 2 + 21a 8 £ 6 ... 



(4) 



We can now make a connection between the above the- 
ory and the band structure. Our calculations, fitted to 
Eq. 2 with n = 3, are presented in Fig. 3. Since the 
high-power coefficients are positive, obviously, renormal- 
ization according to Eq. 3 will lead to a reduction of the 
magnetic susceptibility, x = 1/52 < l/»2- The magni- 
tude of this effect depends on the r.m.s. amplitude, £, of 
the spin fluctuations, which in turn depends on how fast 
x(q) changes at small q's. 

In order to find the value of £ necessary to renormalizc 
the zero-field value of x, one can use Eq. 2 with the n < 
3 expansion: 

dM _ 5 „ 2 35 M 

x ~ ( } = ~dH " fl2 + + y a6? - (5) 

The fit coefficients are a 2 = 478 T/fi B , a 4 = 8990 T/fj,%, 
and a 6 = 277 T//i 5 B . Setting x(0) equal to the experi- 
mental value [30,31] leads to £ — 0.15/ib- However, it is 
highly desirable to find a way of estimating £ in a real 
material using ab initio calculations. This can be done 



using the fluctuation-dissipation theorem along the lines 
suggested by Moriya [33] and elaborated by many authors 
(see, e.g., Refs. [34-36]), which states that for zero-point 
fluctuations 



dui 1 
2^2 



Imx(q,w), 



(6) 



where f2 is the Brillouin zone volume [37] . It is customary 
to approximate x(q> w ) near a QCP as 



X 



^(q, lu) = x \0, 0)-I + cq 2 - iw/Tq, 



(7) 



where Xo (0)0) = ^-/N(E F ) (density of states per spin) 
is the bare (noninteracting) static uniform susceptibility, 
and I is the Stoner parameter which is weakly dependent 
on q and w. Obviously, Xo 1 ( t l' w ) = Xo^O, 0) + c 1 2 ~ 
iui/Tq is the noninteracting susceptibility. Although not 
necessary [35], a convenient approximation, good near a 
QCP, is that x _1 (0,0) w 0, that is, I « l/N(E F ). One 
can also use an expansion for xo(q, w), equivalent to Eq. 
7, namely 



Xo(q, w) = N(E F ) - aq 2 + ibto/q. 



(8) 



Moriya mentioned in his book [33] that the coefficients 
a and b are related, in some approximation, to the band 
structure, in particular, to the effective mass of electrons 
at the Fermi level and to some contour integral along 
a line on the Fermi surface. While Moria's expressions 
are difficult to evaluate numerically within the standard 
band structure calculations, one can rewrite equivalent 
expressions, better suited for actual calculations. For 
completness, we present below the full derivation: 



Re X o(q,0) - ]T [f(E k ) - /(£ k+q )] (£ k+q - E^ 1 

k 

Im X0 (q,w) = ]T[/(£k) - /(£ k+q )]<K£ k+q -Ek- w), 



(9) 
(10) 



where f(E) is the Fermi function, -^jf 2 = S(E-E F ). Expanding Eq. 9 in A = E^ +Cl -E^ = v k -q+i J2 a f3 M k **9/3 + 

a/3 



we get to second order in q 
Rcxo(q,0)=7V(£ F )+^ 



lfdS{e k -E F )\ E ,/3Mk 9*9/3 1 fd 2 S(e k -E F ) . 
2 { dE F ) (^+— 2 } + 6 { dE 2 F ) (Vk ' q) 



The odd powers of v k cancel out and we get (a, (3 = x,y, z) 

Rpv y N( F \ I ^q a qp d{N(E F )^) ^q aqp d 2 (N(E F )v a v p ) 
Rexo(q) - N(E F ) + ^ — ^ + ^ ~ ^fS 



a, 



N(Ep) + q 2 d(N(E F )f, xx ) + q 2 d 2 (N(E F )v 2 x ) 



dE F 



6 dE 2 F 



(11) 

(12) 
(13) 
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where v 2 = v 2 = v z , ji xx — n yy — [i zz . The last equal- 
ity assumes cubic symmetry; generalization to a lower 
symmetry is trivial. Using the following relation, 



dF(s k ) 



dF (gk) 



¥ VkF,Sk) =¥ *> 

one can prove that 

d 2 (N(E F )v 2 x ) d{N(E F )^ xx ) 



de-L 



Therefore 

Re Xo {cfi=N(E F ) 



^ <p(N(E F )vZ) 
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Similarly for Eq. 10 one has 

dm 



Imxo(q,w) = 



w<5(vk-q — w) 



(14) 



(15) 



(16) 



After averaging over the directions of q, this becomes, 
for small u>, 



Z v^q zq 



vl+vl + v 2 z . 



(17) 



Although the Fermi velocity is obviously different along 
different directions, it is still a reasonable approximation 
to introduce an average v F . Then the frequency cutoff in 
Eq. 17 is uj c « qv F . 

From Eq. 8 it follows that 



Imx(q,w) = 



bquN(E F ) 2 



a 2 q 6 + 6 2 w 2 ' 
and, performing the integrations, 



(18) 



, 2 = bv 2 F N{E F f 
2a 2 n 

_ 3b(N(E F )v 2 x )N(E F ) 



2a 2 n 



[Q 4 ln(l + Q- 4 )+ln(l + Q 4 )] 

[Q 4 ln(l + Q- 4 )+ln(l + Q 4 )], 



where Q = q c 



with q c the cutoff in the momen- 



tum space. There is no solid prescription to estimate 
the cutoff value. At small Q the dependence of £ on Q is 
quadratic, however, at large Q it becomes relatively weak 
(logarithmic). While the susceptibility X (<i,oj) can > m 
principle, be calculated exactly, there is no rigorous def- 
inition of q c . The conceptual difficulty here is, as in all 
problems related to electron-electron interactions, that 
some part of the effect in question is already included in 
the LDA, and rigorous treatment of the double-counting 
becomes virtually impossible (cf. discussion of this issue 
in connection to the LDA+U method [39]). At this point 



one needs to make some choice of q c . A natural ansatz is 
to choose the value of q at which the model susceptibility 

(Eq. 15) becomes unphysical (negative), q c = \J ' n (Ef) _ 

The above formulas reduce all parameters needed for 
estimating the r.m.s. amplitude of spin fluctuations 
to four integrals over the Fermi surface: N(E F ), a — 

12 dEj b - 2\ ly \^F) V /' V F - V N(E F ) ' 

It should be noted that these integrals are extremely sen- 
sitive to the k-point mesh. We used various meshes be- 
tween 40x40x40 and 60x60x60, and averaged the results 
using the bootstrap method [41] (to eliminate the effect 
of special points coinciding with mesh points). Velocities 
were calculated [42] as matrix elements of the momentum 
operator, using the optic program of the WIEN package. 
We obtained (all energies are measured in Ry, lengths 
in Bohr, and velocities in Ry-Bohr) N(E F ) = 17.1, 

(N{E F )v 2 x ) = 0.58, d2 W| a *>'> = 1700, (N{E F )v^) 



(N(E F )vl) _ 
N(E F ) 



= 0.31. Correspondingly, a ■■ 



= 135, v F - 

140, b w 72, and q c = \J N ^ F ' > = 0.35, using the above- 
mentioned ansatz. 
Now we get 

£ = 0.2 MB v/Q 4 ln(l + Q- 4 ) + ln(l + Q 4 ), (19) 

and with Q = 0.88, we obtain £ =0.16 \lb- Note that 
the energy of a long-range spin fluctuation with such an 
amplitude is of the order of a few meV per atom, as can 
be seen from Fig. 1. 

This result is quite sensitive to the second derivative 

d 2 (N(E F )vl) 

— dEi — , which was the most difficult quantity to 

calculate. An inspection of the energy dependence of 
(N{E F )v 2 ) (Fig. 6, inset) elucidates the reason: the 
Fermi energy in Pd lies near an inflection point. As a 

., d 2 (N(E F )vl) . „/-,,, , , ,. 

result, — dE i is small (and hard to calculate reli- 
ably). This, perhaps, is not accidental; were this deriva- 
tive 2-3 times larger, the mean amplitude of spin fluctua- 
tion would have been relatively small even given extreme 
proximity of this material to the ferromagnetic instabil- 
ity, because the relevant phase space would have been 
too small. If this approximation is correct, this gives 
an important hint for identifying quantum critical mate- 
rials from the LDA calculations: the calculated ground 
state should be close to ferromagnetic instability (on ei- 
ther side) and the Fermi energy should be close to an 
inflection point of the (N(E)v 2 ,). 

The calculated value of £, if substituted into Eq. 3, 
gives x ~ 6-4 x 10~ 4 emu/mol, practically the same as 
the experimenatal number. Such a good agreement is 
without doubt fortuitous; for instance, using GGA as a 
starting point instead of LDA would have destroyed this 
agreement. [40] We should keep in mind that, first of all, 
the formalism itself is very crude; Xo (q, was expanded 
to leading terms at small q, but this expansion is used up 
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to some large q c comparable with hp. Furthermore, a key 
parameter in the formalism is the cut-off momentum q c , 
for which we use an ansatz based on the large-q behavior 
of the model x(g, uf). 

However, the fact that this procedure produces a cor- 
rection of the right order of magnitude is probably robust 
and suggests that the underlying physics was identified 
correctly. 

To summarize, we use highly accurate LDA calcula- 
tions to estimate the parameters in Moriya's spin fluctu- 
ation theory, and thereby estimate the corrections, due 
to long wavelength spin fluctuations, to the LDA results. 
Let us, in conclusion, repeat our main points. The key 
parameter defining the nontrivial physics near the QCP 
is the mean-square amplitude of the spin fluctuations. 
This parameter is a highly material dependent, nonlocal 
quantity, determined by the spin susceptibility in a large 
part of the Brillouin zone, as well as by the characteristic 
cut-off length separating "non-trivial" spin fluctuations 
from spin-fluctuation implicitly included in the LDA. It 
is hoped, however, that this parameter is mainly defined 
by the long wavelength part of the susceptibility, while 
the short wave-length characteristics, including the cut- 
off length, may be only weakly material, pressure, etc., 
dependent. We implements this idea, relating, in the cor- 
responding approximation, the mean-square amplitude 
of the spin fluctuations near a QCP with characteris- 
tics of the one-electron band structure. The formalism 
is based on the (1) Stoncr theory for spin susceptibility, 
(2) fluctuation-dissipation theorem, and (3) lowest-order 
expansion of the real and imaginary part of the polar- 
ization operator in terms of the frequency and the wave 
vector. The actual band structure of the material is taken 
into account via the lowest-order expansion coefficients 
of the LDA susceptibility, while the effects beyond the 
lowest order in q and u arc neglected. Together with 
the Landau expansion of the free energy, also computable 
within the LDA formalism, this allows one to treat quan- 
tum criticality semi-quantitatively on the basis of LDA 
calculations. 
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FIG. 1. Calculated LSDA total energy, E, (in eV) with 
respect to M = 0/jb as a function of calculated magnetic mo- 
ment, M (in hb)- 
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FIG. 2. Applied external magnetic field, H, (in Tesla) as 
a function of the calculated LSDA magnetic moment, M (in 
/is). The total moment is shown together with spin compo- 
nent and the orbital component. 
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FIG. 3. The external magnetic field, H, (in Tesla) as a func- 
tion of the calculated magnetic moments, M (in hb)- The fit 
is to n < 3 in Eqn. 2 
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FIG. 4. Magnetic susceptibility, \, (in states/eV-cell) cal- 
culated from the fit of H, \ = (ff?) _1 > snown as a function 
of M. The dashed line at 21 states/eV-cell corresponds to the 
experimental value of x for Pd [30,31]. 




FIG. 5. Magnetic susceptibility, \, (in states/eV-cell) cal- 
culated from the fit of H, \ = (fl?) -1 ' snown as a function 
of H. The dashed line at 21 states/eV-cell corresponds to the 
experimental value of x for Pd [30,31]. 
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FIG. 6. d ( N< ^' >v *) as a function of energy, calculated by 
the bootstrap method. Note the numerical noise of up to 10%. 
Inset: (N{E F )vl). 
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